** Clearing Stata memory
capture log close
clear all
set more off, perm
set seed 1234

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////// Table O.49: Quantile Regressions: Log (Average Annual Wages) Seven to 12 Years After Admission Exam //////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

** Opening Phase 2 norm_scores dataset 
use "Work Data/Gender_Phase2_long.dta",clear

*** Creating variables
encode subject, gen (sub)
tab subject, gen (d_sub)
label var sub "Subject"

** Subject dummies
rename d_sub1 Biology
rename d_sub2 Chemistry
rename d_sub3 Geography
rename d_sub4 History
rename d_sub5 Language
rename d_sub6 Mathematics
rename d_sub7 Physics
rename d_sub8 Portuguese
* Labels
label var Biology "Biology"
label var Chemistry "Chemistry"
label var Geography "Geography"
label var History "History"
label var Math "Mathematics"
label var Physics "Physics"
label var Portuguese "Portuguese"
label var Language "Foreign Language"

** Interaction: priority X female
gen fem_priority=female*priority
label var fem_priority "Female $\times$ Priority"

** Interaction: priority X subject
foreach v of varlist Biology-Portuguese {
gen fem_`v'=`v'*female
label var fem_`v' "Female $\times$ `v'"
gen prio_`v'=priority*`v'
label var prio_`v' "Priority $\times$ `v'"
gen fem_prio_`v'=fem_priority*`v'
label var fem_prio_`v' "Female $\times$ Priority $\times$ `v'"
}

global subject "Chemistry Geography History Mathematics Physics"
global subject_fem "fem_Chemistry fem_Geography fem_History fem_Mathematics fem_Physics"

** P1 scores: P1 normalized subject-specific scores
forvalues i=2(1)4 {
gen norm_p1score`i'=norm_p1score^`i'
sum norm_p1score`i'
}

*********************************************************************************
****************   Relative performances ****************************************
*********************************************************************************

******************************** ENEM ********************************

foreach v in norm_enem_w {
bys year female: egen `v'_ave_g=mean(`v')
gen `v'_g=`v'-`v'_ave_g
bys year female: sum `v'_g
}
drop norm_enem_w_ave_g

* Priority x relative performance in ENEM:
foreach v in norm_enem_w {
gen `v'_priority_g=`v'_g*priority
forvalues i=2(1)4 {
gen `v'_priority_g`i'=`v'_g^`i'*priority
sum `v'_priority_g`i'
}
}

global g_norm_enem_w_prio norm_enem_w_priority_g*
d $g_norm_enem_w_prio

* Interaction: subject X ENEM
foreach v of varlist Biology-Portuguese {
gen enem_`v'=`v'*norm_enem_w_g
label var enem_`v' "ENEM $\times$ `v'"
gen fem_enem_`v'=female*norm_enem_w_g*`v'
label var fem_enem_`v' "Female $\times$ ENEM $\times$ `v'"
forvalues i=2(1)4 {
gen enem_`v'_`i'=enem_`v'^`i'
gen fem_enem_`v'_`i'=fem_enem_`v'^`i'
}
sum enem_`v'* fem_enem_`v'*
}

global g_pol_enem_sub "enem_Chemistry* enem_Geography* enem_History* enem_Mathematics* enem_Physics*"
d $g_pol_enem_sub

******************************** Phase 1 scores ********************************

foreach v in norm_p1score {

tab year, sum(`v')
bys year subject female: egen gs_`v'_ave=mean(`v')
gen gs_`v'=`v'-gs_`v'_ave
bys year female subject: sum gs_`v'
drop gs_`v'_ave

forvalues i=2(1)4 {
gen gs_`v'`i'=gs_`v'^`i'
sum gs_`v'`i'
}

global gs_pol_`v' gs_`v' gs_`v'2 gs_`v'3 gs_`v'4
d $gs_pol_`v'

* Priority x Phase 1 scores:
gen gs_`v'_prio=gs_`v'*priority
forvalues i=2(1)4 {
gen gs_`v'_prio`i'=gs_`v'`i'*priority
sum gs_`v'_prio*
}
}

global gs_pol_norm_p1score_prio gs_norm_p1score_prio*
d $gs_pol_norm_p1score_prio

*********************************************************************************
**************** Main sample ****************************************************
*********************************************************************************

* 1) Only years before the affirmative action took place
drop if aa_year==1
tab year
drop if year==2000
tab year

* 2) Drop Portuguese and Foreign Language (in Phase 1 there is no Portuguese or Foreign Language exams - For Portuguese Phase 1 has an essay)
 tab subject, sum(norm_p1score)
 drop if subject=="lang" | subject=="port" 
 tab subject, sum(norm_p1score)
 drop Language Portuguese prio_Language prio_Portuguese fem_prio_Language fem_prio_Portuguese 

 ******************************** SAVE DATA SET TO BOOTSTRAP ********************************
 
keep norm_score fem_priority priority $subject $subject_fem $gs_pol_norm_p1score $g_norm_enem_w_prio $gs_pol_norm_p1score_prio inscri2 female ///
enem norm_enem_w gen_ques_st1 essay_st1 total_st1 career_choice year $g_pol_enem_sub

save "Work Data/Temp_Bootstrap.dta", replace

********************************************************************************************************
****************  Step 1: In the first step we obtain initial estimates and store them  ****************
********************************************************************************************************

** # 1 - Main specification, exclude interaction term
reghdfe norm_score priority $subject $subject_fem $g_pol_enem_sub $g_norm_enem_w_prio $gs_pol_norm_p1score $gs_pol_norm_p1score_prio, cluster(inscri2) absorb(inscri2) resid
** # 2 - Save residuals 
predict residuals, resid
** # 3 - Create measure of relative priority performance
ttest residuals, by(female) // Residuals close to zero in both cases
ttest residuals if priority==1, by(female)
sum residuals if female==1 & priority==1
bys inscri2: egen res_prio_ave=mean(residuals) if priority==1
bys inscri2: egen res_prio_average=min(res_prio_ave)
mdesc res_prio_average
tab career_choice if  res_prio_average==.
bys inscri2: egen res_nonprio_ave=mean(residuals) if priority==0
bys inscri2: egen res_nonprio_average=min(res_nonprio_ave)
mdesc res_nonprio_average
gen res_diff=res_prio_average-res_nonprio_average
mdesc res_diff
ttest res_diff, by(female)

collapse (mean) res_diff female enem norm_enem_w gen_ques_st1 essay_st1 total_st1 career_choice /* expected_duration* */ year, by(inscri2)

count if res_diff==.

 ***************************************************
 ****************** Wages RAIS *********************
 ***************************************************

count
merge 1:1 inscri2 using "Work Data\RAIS_cleaned.dta"
drop if _merge==2
tab _merge  

** # 4 - Wage measures: 

forvalues i=7(1)12 {
gen yearafter`i'=year+`i'
tab yearafter`i', mi
gen annual_wage_after`i'=.
levelsof yearafter`i', local(levels) 
foreach l of local levels {
replace annual_wage_after`i'=mwagetot`l' if yearafter`i'==`l'
}
}

**** Average wage

foreach x in annual_wage {
    
* Average
egen avg_`x'_712=rowmean(`x'_after7 `x'_after8 `x'_after9 `x'_after10 `x'_after11 `x'_after12) // 7-12 years after admission exam
sum avg_`x'*
sum avg_`x'* if _merge==1

}

* Log variables
foreach v of varlist avg*  {
 gen l_`v'=log(`v')
 sum l_`v'
}

foreach v of varlist res_diff {
sum `v'
gen mean_`v'=r(mean)
gen sd_`v'=r(sd)
gen norm_`v'=(`v'-mean_`v')/sd_`v'
sum  norm_`v'
}


******************************** SAVE WAGES DATA TO MERGE WITH DATA SET TO BOOTSTRAP ********************************
preserve
keep l_*_annual_wage_712 inscri2 res_diff norm_res_diff career_choice norm_enem_w female year
renvarlab res_diff norm_res_diff career_choice norm_enem_w female, postfix(_orig)  
duplicates report
save "Work Data/Temp_Wages_Bootstrap.dta", replace
restore

** # 5 - Control for the residuals 
drop if res_diff==.

**********  Annual Wages  **********

foreach v of varlist l_avg_annual_wage_712 {

foreach q in 10 25 50 75 90 {

preserve

capture erase "Work Data\Bootstrap_`v'_q`q'.dta"

** Do not control for program FE

rifhdreg `v'  i.female, over(female) rif(q(`q'))  vce(robust) absorb(year)
gen b_female1=_b[1.female]
gen se_female1=_se[1.female]
foreach x in _cons {
gen b_`x'1=_b[`x']
gen se_`x'1=_se[`x']	
}

rifhdreg `v' i.female norm_res_diff, over(female) rif(q(`q'))  vce(robust) absorb(year)
gen b_female2=_b[1.female]
gen se_female2=_se[1.female]
foreach x in _cons norm_res_diff {
gen b_`x'2=_b[`x']
gen se_`x'2=_se[`x']	
}

rifhdreg `v' i.female norm_enem_w, over(female) rif(q(`q'))  vce(robust) absorb(year)
gen b_female3=_b[1.female]
gen se_female3=_se[1.female]
foreach x in _cons norm_enem_w  {
gen b_`x'3=_b[`x']
gen se_`x'3=_se[`x']	
}

rifhdreg `v' i.female norm_res_diff norm_enem_w,  over(female) rif(q(`q')) vce(robust) absorb(year)
gen b_female4=_b[1.female]
gen se_female4=_se[1.female]
foreach x in _cons norm_res_diff norm_enem_w  {
gen b_`x'4=_b[`x']
gen se_`x'4=_se[`x']	
}

** Control for program FE

rifhdreg `v' i.female i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female5=_b[1.female]
gen se_female5=_se[1.female]
foreach x in _cons  {
gen b_`x'5=_b[`x']
gen se_`x'5=_se[`x']	
}

rifhdreg `v' i.female norm_res_diff i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female6=_b[1.female]
gen se_female6=_se[1.female]
foreach x in _cons  norm_res_diff {
gen b_`x'6=_b[`x']
gen se_`x'6=_se[`x']	
}

rifhdreg `v' i.female norm_enem_w i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female7=_b[1.female]
gen se_female7=_se[1.female]
foreach x in _cons norm_enem_w  {
gen b_`x'7=_b[`x']
gen se_`x'7=_se[`x']	
}

rifhdreg `v' i.female norm_enem_w norm_res_diff i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female8=_b[1.female]
gen se_female8=_se[1.female]
foreach x in _cons norm_res_diff norm_enem_w  {
gen b_`x'8=_b[`x']
gen se_`x'8=_se[`x']	
}

gen boot_run=0
keep boot_run b_* se_*
keep if _n==1  
count
sum

save "Work Data/Bootstrap_`v'_q`q'.dta", replace

restore

}
}

********************************************************************************************************
**************************************  Step 2: Bootstrap   ********************************************
********************************************************************************************************

use "Work Data/Temp_Bootstrap.dta", clear
merge m:1 inscri2 using "Work Data/Temp_Wages_Bootstrap.dta"
drop _merge
count 

local boots = 999

quietly {
forvalues i=1(1)`boots' {

noisily display "Working on `i' out of `boots' at $S_TIME"

foreach v of varlist l_avg_annual_wage_712 {

foreach q in 10 25 50 75 90 {

preserve

local seed=2023 + `i'
set seed `seed'

bsample,  cluster(inscri2) idcluster(newid)  

********** Estimating the residuals **********

** # 1 - Main specification, exclude interaction term
reghdfe norm_score priority $subject $subject_fem $g_pol_enem_sub $g_norm_enem_w_prio $gs_pol_norm_p1score $gs_pol_norm_p1score_prio, cluster(newid) absorb(newid) resid
** # 2 - Save residuals 
predict residuals, resid
** # 3 - Create measure of relative priority performance
bys newid: egen res_prio_ave=mean(residuals) if priority==1
bys newid: egen res_prio_average=min(res_prio_ave)
mdesc res_prio_average
bys newid: egen res_nonprio_ave=mean(residuals) if priority==0
bys newid: egen res_nonprio_average=min(res_nonprio_ave)
mdesc res_nonprio_average
gen res_diff=res_prio_average-res_nonprio_average
mdesc res_diff
ttest res_diff, by(female) 

collapse (mean) res_diff female enem norm_enem_w gen_ques_st1 essay_st1 total_st1 year career_choice `v', by(newid)

drop if res_diff==. 
count if res_diff==.
ttest res_diff, by(female)

* Normalizing residuals
foreach r of varlist res_diff {
sum `r'
gen mean_`r'=r(mean)
gen sd_`r'=r(sd)
gen norm_`r'=(`r'-mean_`r')/sd_`r'
sum  norm_`r'
}

** Do not control for program FE

rifhdreg `v'  i.female, over(female) rif(q(`q'))  vce(robust) absorb(year)
gen b_female1=_b[1.female]
foreach x in _cons {
gen b_`x'1=_b[`x']
}

rifhdreg `v' i.female norm_res_diff, over(female) rif(q(`q'))  vce(robust) absorb(year)
gen b_female2=_b[1.female]
foreach x in _cons norm_res_diff {
gen b_`x'2=_b[`x']
}

rifhdreg `v' i.female norm_enem_w, over(female) rif(q(`q'))  vce(robust) absorb(year)
gen b_female3=_b[1.female]
foreach x in _cons norm_enem_w  {
gen b_`x'3=_b[`x']
}

rifhdreg `v' i.female norm_res_diff norm_enem_w,  over(female) rif(q(`q')) vce(robust) absorb(year)
gen b_female4=_b[1.female]
foreach x in _cons norm_res_diff norm_enem_w  {
gen b_`x'4=_b[`x']
}

** Control for program FE

rifhdreg `v' i.female i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female5=_b[1.female]
foreach x in _cons  {
gen b_`x'5=_b[`x']
}

rifhdreg `v' i.female norm_res_diff i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female6=_b[1.female]
foreach x in _cons  norm_res_diff {
gen b_`x'6=_b[`x']
}

rifhdreg `v' i.female norm_enem_w i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female7=_b[1.female]
foreach x in _cons norm_enem_w  {
gen b_`x'7=_b[`x']
}

rifhdreg `v' i.female norm_enem_w norm_res_diff i.year, over(female) rif(q(`q'))  vce(robust) absorb(career_choice) 
gen b_female8=_b[1.female]
foreach x in _cons norm_res_diff norm_enem_w  {
gen b_`x'8=_b[`x']
}

gen boot_run=`i'
keep boot_run b_*
keep if _n==1  
summarize

append using "Work Data/Bootstrap_`v'_q`q'.dta" 
sleep 50
save "Work Data/Bootstrap_`v'_q`q'.dta", replace

restore
}

}
}
}

********************************************************************************************************
********************************************************************************************************
********************************************************************************************************

estimates clear

foreach q in 10 25 50 75 90 {
	
foreach y of varlist l_avg_annual_wage_712 {
	
********************************************************************************************************
*********************** Step 3: Generating the bootstrapped standard errors ****************************
********************************************************************************************************

preserve	
	
use "Work Data/Bootstrap_`y'_q`q'.dta", clear

** Calculating bootstrapped standard errors

foreach x in female norm_enem_w norm_res_diff _cons {
forvalues i = 1(1)8 {
capture egen BS_sd_`x'`i'=sd(b_`x'`i')  if boot_run~=0
}
}

** Making a matrix with bootstrapped standard errors

matrix define bsemat`y'1 = J(1,3,0)
matrix define bsemat`y'2 = J(1,4,0)
matrix define bsemat`y'3 = J(1,4,0)
matrix define bsemat`y'4 = J(1,5,0)

matrix colnames bsemat`y'1 = 0b.female 1.female _cons 
matrix colnames bsemat`y'2 = 0b.female 1.female resid _cons
matrix colnames bsemat`y'3 = 0b.female 1.female norm_enem_w _cons
matrix colnames bsemat`y'4 = 0b.female 1.female resid norm_enem_w _cons

matrix define bsemat`y'5 = J(1,7,0)
matrix define bsemat`y'6 = J(1,8,0)
matrix define bsemat`y'7 = J(1,8,0)
matrix define bsemat`y'8 = J(1,9,0)

matrix colnames bsemat`y'5 = 0b.female 1.female 2001b.year 2002.year 2003.year 2004.year _cons
matrix colnames bsemat`y'6 = 0b.female 1.female resid 2001b.year 2002.year 2003.year 2004.year _cons
matrix colnames bsemat`y'7 = 0b.female 1.female norm_enem_w 2001b.year 2002.year 2003.year 2004.year _cons
matrix colnames bsemat`y'8 = 0b.female 1.female resid norm_enem_w 2001b.year 2002.year 2003.year 2004.year _cons
			  
forvalues i = 1(1)8 {
sum BS_sd_female`i'
scalar BS_sd_female = r(mean)
matrix bsemat`y'`i'[1,2] = BS_sd_female
}

foreach i in 2 4 6 8 {
sum BS_sd_norm_res_diff`i'
scalar BS_sd_norm_res_diff = r(mean)
matrix bsemat`y'`i'[1,3] = BS_sd_norm_res_diff
}


foreach i in 3 7 {
sum BS_sd_norm_enem_w`i'
scalar BS_sd_norm_enem_w = r(mean)
matrix bsemat`y'`i'[1,3] = BS_sd_norm_enem_w
}

foreach i in 4 8 {
sum BS_sd_norm_enem_w`i'
scalar BS_sd_norm_enem_w = r(mean)
matrix bsemat`y'`i'[1,4] = BS_sd_norm_enem_w
}

foreach i in 1 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,3] = BS_sd_cons
}

foreach i in 2 3 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,4] = BS_sd_cons
}

foreach i in 4 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,5] = BS_sd_cons
}


foreach i in 5 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,7] = BS_sd_cons
}

foreach i in 6 7 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,8] = BS_sd_cons
}

foreach i in 8 {
sum BS_sd__cons`i'
scalar BS_sd_cons = r(mean)
matrix bsemat`y'`i'[1,9] = BS_sd_cons
}

forvalues i = 1(1)8 {
matrix list bsemat`y'`i'
}

restore

********************************************************************************************************
*********************** Step 4: Computing p-values and exporting tables ********************************
********************************************************************************************************

use "Work Data/Temp_Wages_Bootstrap.dta", clear
rename norm_res_diff_orig resid
rename career_choice_orig career_choice
rename norm_enem_w_orig norm_enem_w
rename female_orig female

label var norm_enem_w "Norm. ENEM scores"
label var female "Female"
label var resid "Relative priority performance"

**********  Annual Wages **********

drop if resid==.

** Do not control for program FE

rifhdreg `y'  i.female, over(female) rif(q(`q'))  vce(robust) absorb(year)
matrix define b`y'1 = e(b)
scalar DF`y'1=e(df_r)
estimates store `y'`q'1
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'1
estadd matrix seboot =  bsemat`y'1

rifhdreg `y' i.female resid, over(female) rif(q(`q')) vce(robust) absorb(year)
matrix define b`y'2 = e(b)
scalar DF`y'2=e(df_r)
estimates store `y'`q'2
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'2
estadd matrix seboot =  bsemat`y'2

rifhdreg `y' i.female norm_enem_w, over(female) rif(q(`q')) vce(robust) absorb(year)
matrix define b`y'3 = e(b)
scalar DF`y'3=e(df_r)
estimates store `y'`q'3
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'3
estadd matrix seboot =  bsemat`y'3

rifhdreg `y' i.female resid norm_enem_w, over(female) rif(q(`q')) vce(robust) absorb(year)
matrix define b`y'4 = e(b)
scalar DF`y'4=e(df_r)
estimates store `y'`q'4
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'4
estadd matrix seboot =  bsemat`y'4

** Control for program FE

rifhdreg `y' i.female i.year, over(female) rif(q(`q')) vce(robust) absorb(career_choice)
matrix define b`y'5  = e(b)
scalar DF`y'5=e(df_r)
estimates store `y'`q'5
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'5
estadd matrix seboot =  bsemat`y'5

rifhdreg `y' i.female resid i.year, over(female) rif(q(`q')) vce(robust) absorb(career_choice)
matrix define b`y'6  = e(b)
scalar DF`y'6=e(df_r)
estimates store `y'`q'6
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'6
estadd matrix seboot =  bsemat`y'6

rifhdreg `y' i.female norm_enem_w i.year, over(female) rif(q(`q')) vce(robust) absorb(career_choice)
matrix define b`y'7 = e(b)
scalar DF`y'7=e(df_r)
estimates store `y'`q'7
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'7
estadd matrix seboot =  bsemat`y'7

rifhdreg `y' i.female resid norm_enem_w i.year, over(female) rif(q(`q')) vce(robust) absorb(career_choice)
matrix define b`y'8 = e(b)
scalar DF`y'8=e(df_r)
estimates store `y'`q'8
sum `y' if e(sample) & female==0
estadd scalar ymean_men = r(mean)
* Adding bootstrapped standard errors to the estimation results.
estimates restore `y'`q'8
estadd matrix seboot =  bsemat`y'8

* Calculating p-values

foreach x in pval pvalue {	
	
matrix define b`x'`y'1 = J(1,3,0)
matrix define b`x'`y'2 = J(1,4,0)
matrix define b`x'`y'3 = J(1,4,0)
matrix define b`x'`y'4 = J(1,5,0)
matrix define b`x'`y'5 = J(1,7,0)
matrix define b`x'`y'6 = J(1,8,0)
matrix define b`x'`y'7 = J(1,8,0)
matrix define b`x'`y'8 = J(1,9,0)

matrix colnames b`x'`y'1 = 0b.female 1.female _cons 
matrix colnames b`x'`y'2 = 0b.female 1.female resid _cons
matrix colnames b`x'`y'3 = 0b.female 1.female norm_enem_w _cons
matrix colnames b`x'`y'4 = 0b.female 1.female resid norm_enem_w _cons
matrix colnames b`x'`y'5 = 0b.female 1.female 2001b.year 2002.year 2003.year 2004.year _cons
matrix colnames b`x'`y'6 = 0b.female 1.female resid 2001b.year 2002.year 2003.year 2004.year _cons
matrix colnames b`x'`y'7 = 0b.female 1.female norm_enem_w 2001b.year 2002.year 2003.year 2004.year _cons
matrix colnames b`x'`y'8 = 0b.female 1.female resid norm_enem_w 2001b.year 2002.year 2003.year 2004.year _cons_cons

}

foreach x in 1 {
forvalues j = 2/3 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

foreach x in 2 3 {
forvalues j = 2/4 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

foreach x in 4 {
forvalues j = 2/5 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

foreach x in 5 {
foreach j in 2 7 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

foreach x in  6 7 {
foreach j  in 2 3 8 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

foreach x in 8 {
foreach j in 2 3 4 9 {
matrix bpval`y'`x'[1,`j']=2*(normal(-abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
matrix bpvalue`y'`x'[1,`j']=(2 * ttail(DF`y'`x', abs(b`y'`x'[1,`j']/bsemat`y'`x'[1,`j'])))
}
}

forvalues i = 1(1)8 {
matrix list b`y'`i'
matrix list bsemat`y'`i'
matrix list bpval`y'`i'
matrix list bpvalue`y'`i'
}

* Adding my p values to the estimation results.
forvalues i=1(1)8 {
estimates restore `y'`q'`i'
estadd matrix pboot = bpval`y'`i'
}

}
}

estadd local program "No": *1 *2  *3 *4
estadd local year "Yes": *1 *2  *3 *4 *5 *6 *7 *8
estadd local program "Yes": *5  *6 *7  *8


*******************************************************************************************************
********************************** Exporting tables ***************************************************
*******************************************************************************************************

use "Work Data/Temp_Wages_Bootstrap.dta", clear

rename norm_res_diff_orig resid
rename norm_enem_w_orig norm_enem_w
rename female_orig female

label var norm_enem_w "Norm. ENEM scores"
label var female "Female"
label var resid "Relative priority performance"

label var norm_enem_w "Norm. ENEM scores"
label var female "Female"
label var resid "Relative priority performance"

label var female "Female"

** AVERAGE

* 7 to 12 years after the admission exam


* Without major FE

esttab l_avg_annual_wage_712101 l_avg_annual_wage_712102 l_avg_annual_wage_712251 l_avg_annual_wage_712252 l_avg_annual_wage_712501 l_avg_annual_wage_712502 l_avg_annual_wage_712751 l_avg_annual_wage_712752 l_avg_annual_wage_712901 l_avg_annual_wage_712902 using "Output\Annual_avgwages_res_diff_ENEM_712_BS_Quantiles.tex" , nogaps booktabs label  collabels(none) se star(* 0.10 ** 0.05 *** 0.01) keep(1.female  resid) replace f cells(b(star pvalue(pboot) fmt(%9.3f)) seboot(par fmt(%9.3f))) refcat(1.female "\\ \multicolumn{11}{l}{\textit{Panel A: Without major FE}} \\", nolabel) stats(year program, fmt(%3s %3s) labels("Exam year FE" "Major FE"))  nomtitle mgroups("10" "25" "50" "75" "90", pattern(1 0 1 0 1 0 1 0 1 0) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span})) order(1.female resid)  substitute("Female=1" "Female")

* With major FE
esttab l_avg_annual_wage_712105 l_avg_annual_wage_712106 l_avg_annual_wage_712255 l_avg_annual_wage_712256 l_avg_annual_wage_712505 l_avg_annual_wage_712506 l_avg_annual_wage_712755 l_avg_annual_wage_712756 l_avg_annual_wage_712905 l_avg_annual_wage_712906 using "Output\Annual_avgwages_res_diff_ENEM_712_BS_Quantiles.tex" , nogaps booktabs label  collabels(none) se star(* 0.10 ** 0.05 *** 0.01) keep(1.female  resid)  f cells(b(star pvalue(pboot) fmt(%9.3f)) seboot(par fmt(%9.3f))) refcat(1.female "\\ \multicolumn{11}{l}{\textit{Panel B: With major FE}} \\", nolabel) stats(year program sep N ymean_men , fmt(%3s %3s %3s  %9.0fc %7.3f ) labels("Exam year FE" "Major FE"  " " "Number of observations" "Mean dependent variable (Men)")) nomtitle append order(1.female resid)  substitute("Female=1" "Female") nonumbers

***** Erase temporary data sets
erase "Work Data/Temp_Bootstrap.dta"
erase "Work Data/Temp_Wages_Bootstrap.dta"





